Tarea 2: Frequentist Inference

CC6104: Statistical Thinking

Integrantes :

Cuerpo Docente:

Fecha límite de entrega:

Índice:

  1. Objetivo
  2. Instrucciones
  3. Referencias
  4. Elaboración de Código

Objetivo

a la segunda tarea del curso Statistical Thinking. Esta tarea tiene como objetivo evaluar los contenidos teóricos de la segunda parte del curso, los cuales se enfocan principalmente en inferencia estadística, diseño de experimentos y test de hipótesis. Si aún no han visto las clases, se recomienda visitar los enlaces de las referencias.

La tarea consta de una parte teórica que busca evaluar conceptos vistos en clases. Seguido por una parte práctica con el fin de introducirlos a la programación en R enfocada en el análisis estadístico de datos.

Instrucciones:

Referencias:

Slides de las clases:

Enlaces a videos de las clases:

Documentación:

Elaboración de Código

En la siguiente sección deberá resolver cada uno de los experimentos computacionales a través de la programación en R. Para esto se le aconseja que cree funciones en R, ya que le facilitará la ejecución de gran parte de lo solicitado.

Para el desarrollo preste mucha atención en los enunciados, ya que se le solicitará la implementación de métodos sin uso de funciones predefinidas. Por otro lado, Las librerías permitidas para desarrollar de la tarea 2 son las siguientes:

# Manipulación de estructuras
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
# Para realizar plots
library(ggplot2)
library(plotly)
## 
## Adjuntando el paquete: 'plotly'
## 
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## 
## The following object is masked from 'package:stats':
## 
##     filter
## 
## The following object is masked from 'package:graphics':
## 
##     layout
# Manipulación de varios plots en una imagen.
library(gridExtra)
## 
## Adjuntando el paquete: 'gridExtra'
## 
## The following object is masked from 'package:dplyr':
## 
##     combine

Elimine eval=FALSE para ejecutar las celdas.

Pregunta 1: Estimadores.

# Function to calculate bias
calculate_bias <- function(estimates, p_true = 0.5) {
  mean(estimates) - p_true
}

# sigma = 0
cte <- rep(0.5, 1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 0)
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_0 <- data.frame(x = x, p_hat = data)
bias_0 <- calculate_bias(data)

# Create individual plots
p1 <- ggplot(df_0, aes(x = x, y = p_hat)) +
  geom_line(color = "blue") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 0)", x = "n", y = "Estimator of p") +
  theme_minimal()

# sigma = 1
cte <- rep(0.5,1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 1)  # sigma = 1
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_1 <- data.frame(x = x, p_hat = data)
bias_1 <- calculate_bias(data)

p2 <- ggplot(df_1, aes(x = x, y = p_hat)) +
  geom_line(color = "red") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 1)", x = "n", y = "Estimator of p") +
  theme_minimal()

# sigma = 2
cte <- rep(0.5,1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 2)  # sigma = 2
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_2 <- data.frame(x = x, p_hat = data)
bias_2 <- calculate_bias(data)

p3 <- ggplot(df_2, aes(x = x, y = p_hat)) +
  geom_line(color = "green") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 2)", x = "n", y = "Estimator of p") +
  theme_minimal()

# sigma = 4
cte <- rep(0.5,1000)
x <- 1:1000
data <- c()

for (n in x) {
  X <- rbinom(n, 1, 0.5)
  epsilon_sigma <- rnorm(1, 0, 4)  # sigma = 4
  p_hat <- epsilon_sigma + mean(X)
  data <- c(data, p_hat)
}
df_4 <- data.frame(x = x, p_hat = data)
bias_4 <- calculate_bias(data)

p4 <- ggplot(df_4, aes(x = x, y = p_hat)) +
  geom_line(color = "purple") +
  geom_hline(yintercept = 0.5, linetype = "dashed", color = "black") +
  labs(title = "Estimator of p (sigma = 4)", x = "n", y = "Estimator of p") +
  theme_minimal()

# Print plots
print(p1)

print(p2)

print(p3)

print(p4)

# Print out the biases for each sigma
cat("Bias for sigma = 0:", bias_0, "\n")
## Bias for sigma = 0: -0.0005238953
cat("Bias for sigma = 1:", bias_1, "\n")
## Bias for sigma = 1: -0.01096729
cat("Bias for sigma = 2:", bias_2, "\n")
## Bias for sigma = 2: 0.01166566
cat("Bias for sigma = 4:", bias_4, "\n")
## Bias for sigma = 4: 0.08243936

Respuesta

Respuesta

Estos muestran alta variabilidad, tanto para altos como bajos valores de \(\sigma\). El único que mantiene cierta regularidad es \(\sigma = 0\). Esto se debe a que para \(\sigma = 0\) no hay ruido en el estimador, por lo que es esperable que este se acerque al valor real a medida que aumenta n. Vemos que los otros valores de sigma si se diferencian del valor real, aumentando esta diferencia promedio a medida que aumenta n (para \(\sigma = 1\) y \(\sigma = 2\) es cercana a \(\pm 4\), mientras que para \(\sigma = 4\) llega hasta \(\pm 15\)). En conclusión, los estimadores con menos sesgo si son más precisos.

Pregunta 2: Intervalo de Confidencia

En las preguntas 2, 3 y 4 deberá trabajar con el dataset diabetes_prediction_dataset.

datos <- read.csv("diabetes_prediction_dataset.csv", header = TRUE)
head(datos)

Una forma de estimar la distribución de la media de una población es a través de la realización de la sampling distribution de una distribución cualquiera. El valor obtenido en cada sampling distribution nos entrega un estimador de la media que posee una determinada distribución de la población. Sabiendo esto, calcule la sampling distribution de las columnas bmi y blood_glucose_level, obteniendo el intervalo de confianza de \(95\%\) para cada una de las medias obtenidas desde la distribución. Para realizar este experimento considere los siguientes puntos:

Hints:

  • Para comparar las 2 estimaciones puede ser útil normalizar bmi y blood_glucose_level a una misma escala.
  • Para realizar la sampling distribution podría serle útil el comando sample().
  • Puede ser útil generar un dataframe para graficar con ggplot2.

Respuesta:

# Definimos tamaños de muestreo
tamano_muestra = 100
n_muestras = 5000

# Inicializamos estructuras necesarias
list_mean_bmi = vector()
list_interval_bmi_left = vector()
list_interval_bmi_right = vector()
sucess_bmi = 0

list_mean_bgl = vector()
list_interval_bgl_left = vector()
list_interval_bgl_right = vector()
sucess_bgl = 0

# Obtenemos la media y desviación estándar de la población real
real_mean_bmi <- mean(datos$bmi)
real_sd_bmi <- sd(datos$bmi)

real_mean_bgl <- mean(datos$blood_glucose_level)
real_sd_bgl <- sd(datos$blood_glucose_level)

# Sampling distribution, calculo del intervalo de confianza y proporción.
for (i in 1:n_muestras) {
  
  # Generating samples
  sample_bmi <- sample(datos$bmi, tamano_muestra)
  sample_bgl <- sample(datos$blood_glucose_level, tamano_muestra)
  
  # Sampling distribution of the mean
  list_mean_bmi <- c(list_mean_bmi, mean(sample_bmi))
  list_mean_bgl <- c(list_mean_bgl, mean(sample_bgl))
  
  # Standard deviations of the samples
  sd_bmi <- sd(sample_bmi)
  sd_bgl <- sd(sample_bgl)
  
  # Calculating bmi confidence interval
  margin_error_bmi <- 1.96 * (sd_bmi / sqrt(tamano_muestra))
  ci_lower_bmi <- mean(sample_bmi) - margin_error_bmi
  ci_upper_bmi <- mean(sample_bmi) + margin_error_bmi
  list_interval_bmi_left <- rbind(list_interval_bmi_left, ci_lower_bmi)
  list_interval_bmi_right <- rbind(list_interval_bmi_right, ci_upper_bmi) 
  
  # Calculating bgl confidence interval
  margin_error_bgl <- 1.96 * (sd_bgl / sqrt(tamano_muestra))
  ci_lower_bgl <- mean(sample_bgl) - margin_error_bgl
  ci_upper_bgl <- mean(sample_bgl) + margin_error_bgl
  list_interval_bgl_left <- rbind(list_interval_bgl_left, ci_lower_bgl)
  list_interval_bgl_right <- rbind(list_interval_bgl_right, ci_upper_bgl)
  
  # Verify if the real mean is within the confidence interval
  if (ci_lower_bmi <= real_mean_bmi && ci_upper_bmi >= real_mean_bmi) {
    sucess_bmi <- sucess_bmi + 1
  }
  
  if (ci_lower_bgl <= real_mean_bgl && ci_upper_bgl >= real_mean_bgl) {
    sucess_bgl <- sucess_bgl + 1
  }
}

# Proporción de intervalos que contienen la media real
prop_bmi <- sucess_bmi / n_muestras
prop_bgl <- sucess_bgl / n_muestras

cat("Proporción de intervalos que contienen la media real (BMI):", prop_bmi, "\n")
## Proporción de intervalos que contienen la media real (BMI): 0.9474
cat("Proporción de intervalos que contienen la media real (Blood Glucose Level):", prop_bgl, "\n")
## Proporción de intervalos que contienen la media real (Blood Glucose Level): 0.953
# Convertir los datos a dataframes
df_bmi <- data.frame(
  mean = list_mean_bmi,
  lower = list_interval_bmi_left,
  upper = list_interval_bmi_right
)

df_bgl <- data.frame(
  mean = list_mean_bgl,
  lower = list_interval_bgl_left,
  upper = list_interval_bgl_right
)

df_bmi <- head(df_bmi, 100)
df_bgl <- head(df_bgl, 100)

# Columna que nos indica si la media real está dentro del intervalo de confianza calculado
df_bmi$outside_ci <- ifelse(df_bmi$lower > real_mean_bmi | df_bmi$upper < real_mean_bmi, "Fuera", "Dentro")
df_bgl$outside_ci <- ifelse(df_bgl$lower > real_mean_bgl | df_bgl$upper < real_mean_bgl, "Fuera", "Dentro")

# Normalización de las columnas
# Esto no lo usé al final pero lo dejé igual
normalize <- function(x) {
  return ((x - min(x)) / (max(x) - min(x)))
}

# Plot de Intervalos de confianza
# Gráfico para BMI
ggplot(df_bmi, aes(y = 1:100, x = mean)) +
  geom_point(aes(color = outside_ci)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper, color = outside_ci), height = 0.2) +
  scale_color_manual(values = c("Dentro" = "black", "Fuera" = "red")) +
  geom_vline(xintercept = real_mean_bmi, color = "red", linetype = "dashed") + 
  labs(title = "Intervalos de Confianza para BMI", x = "Valor", y = "Muestra") +
  theme_minimal()

# Plot de Intervalos de confianza
# Gráfico para Blood Glucose Level
ggplot(df_bgl, aes(y = 1:100, x = mean)) +
  geom_point(aes(color = outside_ci)) +
  geom_errorbarh(aes(xmin = lower, xmax = upper, color = outside_ci), height = 0.2) +
  scale_color_manual(values = c("Dentro" = "black", "Fuera" = "red")) +
  geom_vline(xintercept = real_mean_bgl, color = "red", linetype = "dashed") +
  labs(title = "Intervalos de Confianza para Blood Glucose Level", x = "Valor", y = "Muestra") +
  theme_minimal()

Respuesta La proporción de intervalos que contienen la media real parece ser muy buena y similar para ambas variables (Cercana al 94%). En general, ambas variables parecen mostrar resultados proporcionalmente similares a pesar de estar en una escala disitinta. Ambas variables suelen mostrar una cantidad baja de intervalos que no contienen la media real, a pesar de estar en una escala de magnitud distinta. Esto se debe a que se están haciendo la misma cantidad de experimentos, y estos al ser muchos, tienden a mostrar una proporción positiva en cuanto a contener la media real dentro del intervalo de confianza.

Pregunta 3: Estimación de Máxima Verosimilitud

El objetivo de esta parte será obtener los parámetros de las distribuciones de la media y la mediana de bmi. Para responder la pregunta realice los siguientes puntos:

Cabe señalar que el método de máxima verosimilitud deberá ser programado por usted y no podrá utilizar librerías que entreguen el valor directo (por ejemplo la librería MASS).

Respuesta

# Obtenemos métricas de muestreo con repetición
means <- c()
medis <- c()
n_muestras <- 10000
tamano_muestra <- 100

for (i in 1:n_muestras) {
  
}

# dataframe con las medias y medianas de las muestras
# Histograma media
# Histograma mediana
# Media
# función log likelihood
ll_plot <- function(mu, sigma) {
  # ...
}

ll_plot <- Vectorize(ll_plot)

# definir espacio donde se va a evaluar ll_plot
mu_m <- 
sigma_m <- 
likelihood_m <- outer(X=mu_m, Y=sigma_m, ll_plot)

filled.contour(x=mu_m, y=sigma_m, z=likelihood_m, xlab=expression(mu), 
               ylab=expression(sigma), color = topo.colors)
# Mediana
# función log likelihood
ll_plot <- function(mu, sigma) {
  # ...
}

ll_plot <- Vectorize(ll_plot)

# definir espacio donde se va a evaluar ll_plot
mu_d <-
sigma_d <-
likelihood_d <- outer(X=mu_d, Y=sigma_d, ll_plot)

filled.contour(x=mu_d, y=sigma_d, z=likelihood_d, xlab=expression(mu), 
               ylab=expression(sigma), color = topo.colors)
likelihood_mean <- function(param) {
  # Definimos los parametros de entrada de la funcion
  mu <- param[1]
  sigma <- param[2]
  # Definimos la likelihood como la suma logaritmica de la función de densidad
  # ...
}

# Agrrgue el rango donde operará nlminb
nlminb(objective=likelihood_mean, start=c(, ), lower=c(, ), upper=c(, ))
likelihood_med <- function(param) {
  # Definimos los parametros de entrada de la funcion
  mu <- param[1]
  sigma <- param[2]
  # Definimos la likelihood como la suma logaritmica de la función de densidad
  # ...
}

# Agrrgue el rango donde operará nlminb
nlminb(objective=likelihood_med, start=c(, ), lower=c(, ), upper=c(, ))
# Muestra de medias


# Histogrmas
# Muestra de medianas


# Histogrmas

Pregunta 4: Bootstrap I

En esta pregunta se error e inetrvalo de confianza para la varianza de la columna HbA1c_level.

Las actividades por realizar son las siguientes:

Respuesta:

# Bootstrap
B <- 5000
largo <- length(datos$HbA1c_level)
output <- c()

for (it in 1:B) {
  
}
# Gráfico de puntos
# Histogrma
# Obtenemos error e intervalos de confianze
se_vars <-
z_a2 <-
alpha <-

# límite inferios
l.CI <- 
# límite superior
u.CI <- 

sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)
sprintf('El SE de la varianzas es: (%.3f)', se_vars)

Respuesta

Pregunta 5: Bootstrap II

El siguiente código genera una regresión lineal de bmi con respecto a age, es decir, una función lineal de age, \(y(age) = b + m\cdot age\), que pretende determinar el valor de bmi.

linearMod <- lm(bmi ~ age, data=datos)
print(linearMod)
## 
## Call:
## lm(formula = bmi ~ age, data = datos)
## 
## Coefficients:
## (Intercept)          age  
##    23.15536      0.09945
m <- linearMod$coefficients["age"]
b <- linearMod$coefficients["(Intercept)"]

ggplot() +
  geom_point(aes(x=datos$age, y=datos$bmi)) +
  geom_line(aes(x=datos$age, y=datos$age*m+b), color="red") +
  ggtitle("Regresión lineal") +
  ylab("bmi") + 
  xlab("Age")

Repita el proceso de la pregunta 4 para estimar el error e intervalos de confianza de los parámetros de la regresión (\(m\) y \(b\)). ¿Qué indican los resultados de Bootstrap?¿Un valor bajo en el error significa que la regresión es buena?

# Bootstrap
B <- 500
largo <- 1000
output_m <- c()
output_b <- c()

for (it in 1:B) {
  
}
# Obtenemos error e intervalos de confianze
se_m <- 
z_a2 <-
alpha <-

# límite inferios
l.CI <- 
# límite superior
u.CI <- 

sprintf('El intervalo de confidencia de 95%% de las varianzas es: (%.3f,%.3f)', l.CI , u.CI)
sprintf('El SE de la varianzas es: (%.3f)', se_m)
# Obtenemos error e intervalos de confianze
se_b <- 
z_a2 <-
alpha <-

# límite inferios
l.CI <- 
# límite superior
u.CI <- 

sprintf('El intervalo de confidencia de 95%% de los sesgos es: (%.3f,%.3f)', l.CI , u.CI)
sprintf('El SE del sesgo es: (%.3f)', se_b)

Respuesta

 

A work by CC6104